!-------------------------------------------------------------LICENSE--------------------------------------------------------------!
!                                                                                                                                  !
!The MAP code is written in Fortran language for magnetohydrodynamics (MHD) calculation with the adaptive mesh refinement (AMR)    !
!and Message Passing Interface (MPI) parallelization.                                                                              !
!                                                                                                                                  !
!Copyright (C) 2012                                                                                                                !
!Ronglin Jiang                                                                                                                     !
!rljiang@ssc.net.cn                                                                                                                !
!585 Guoshoujing Road. Pudong, Shanghai, P.R.C. 201203                                                                             !
!                                                                                                                                  !
!This program is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License         !
!as published by the Free Software Foundation; either version 2 of the License, or (at your option) any later version.             !
!                                                                                                                                  !
!This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of    !
!MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more details.                        !
!                                                                                                                                  !
!You should have received a copy of the GNU General Public License along with this program; if not, write to the Free Software     !
!Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.                                                   !
!                                                                                                                                  !
!-------------------------------------------------------------LICENSE--------------------------------------------------------------!

!==================================================================================================================================|
subroutine bnd_user(ro, mx, my, mz, bx, by, bz, en, gx, gy, gz, po,                                                                &
                    ro0, mx0, my0, mz0, bx0, by0, bz0, en0,                                                                        &
                    x, y, z, nx, ny, nz, dir, sides)
!==================================================================================================================================|
!  sides = 'l' for lower boundary
!  sides = 'h' for higher boundary
!  sides = 'b' for both
!----------------------------------------------------------------------------------------------------------------------------------|

   use parameters
   implicit none

   character(1), intent(in) ::  dir, sides
   integer(4), intent(in) :: nx, ny, nz
   real(8), intent(in) :: x(nx), y(ny), z(nz)
   real(8), dimension(nx, ny, nz), intent(inout) :: ro, mx, my, mz, bx, by, bz, en, gx, gy, gz, po
   real(8), dimension(nx, ny, nz), intent(in) :: ro0, mx0, my0, mz0, bx0, by0, bz0, en0

   integer(4) :: i_s, j_s, k_s, i_e, j_e, k_e
!   real(8) :: sin_theta, cos_theta, sin_phi, cos_phi, a_phi, w, x0, y0, z0, q, a, bt, r0, r, r_xy, br, btheta, bphi, y_p
   real(8) :: t_emerge_end, pr, r0, bt, e_dir, r, x0, y0, z0

!----------------------------------------------------------------------------------------------------------------------------------|
!  array index
!----------------------------------------------------------------------------------------------------------------------------------|
   if (dir .eq. 'x' .and. sides .eq. 'l') then
      i_s = 1
      i_e = ng
      j_s = 1
      j_e = ny
      k_s = 1
      k_e = nz
   elseif (dir .eq. 'x' .and. sides .eq. 'h') then
      i_s = nx - ng + 1
      i_e = nx
      j_s = 1
      j_e = ny
      k_s = 1
      k_e = nz
   elseif (dir .eq. 'y' .and. sides .eq. 'l') then
      i_s = 1
      i_e = nx
      j_s = 1
      j_e = ng
      k_s = 1
      k_e = nz
   elseif (dir .eq. 'y' .and. sides .eq. 'h') then
      i_s = 1
      i_e = nx
      j_s = ny - ng + 1
      j_e = ny
      k_s = 1
      k_e = nz
   elseif (dir .eq. 'z' .and. sides .eq. 'l') then
      i_s = 1
      i_e = nx
      j_s = 1
      j_e = ny
      k_s = 1
      k_e = ng
   elseif (dir .eq. 'z' .and. sides .eq. 'h') then
      i_s = 1
      i_e = nx
      j_s = 1
      j_e = ny
      k_s = nz - ng + 1
      k_e = nz
   endif

!----------------------------------------------------------------------------------------------------------------------------------|
!  emerging flux
!----------------------------------------------------------------------------------------------------------------------------------|
   x0 = 0.0d0
   y0 = 0.0d0
   z0 = 0.0d0
   t_emerge_end = 20.0d0
!   y_p = 40.0d0

   if (t .le. t_emerge_end) then
      y0 = -0.0d0 !- y_p * (1.0d0 - t / t_emerge_end)
      bt = 5.0d0 !* t / t_emerge_end
!      v_e = r / t_emerge_end
   else
      y0 = -0.0d0
      bt = 5.0d0
!      v_e = r / t_emerge_end * (2.0d0 - t / t_emerge_end)
   endif

   r0 = 8.0d0
   e_dir = 1

!   q = -0.1d0             ! twisting number
!   a = 0.04d0 * 200
!   r0 = 0.08d0 * 50

!   do k = k_s, k_e
!   do j = j_s, j_e
!   do i = i_s, i_e
   do k = 1, nz
   do j = 1, ny
   do i = 1, nx

      ro(i, j, k) = ro0(i, j, k)
      mx(i, j, k) = mx0(i, j, k)
      my(i, j, k) = my0(i, j, k)
      mz(i, j, k) = mz0(i, j, k)
      bx(i, j, k) = bx0(i, j, k)
      by(i, j, k) = by0(i, j, k)
      bz(i, j, k) = bz0(i, j, k)
      en(i, j, k) = en0(i, j, k)
      po(i, j, k) = 0.0d0

      pr = (en(i, j, k) -                                                                                                          &
           (mx(i, j, k) * mx(i, j, k) + my(i, j, k) * my(i, j, k) + mz(i, j, k) * mz(i, j, k)) / ro(i, j, k) / 2.0d0 -             &
           (bx(i, j, k) * bx(i, j, k) + by(i, j, k) * by(i, j, k) + bz(i, j, k) * bz(i, j, k)) / 2.0d0) * gmm1

      r = sqrt((x(i) - x0) * (x(i) - x0) + (y(j) - y0) * (y(j) - y0) + (z(k) - z0) * (z(k) - z0))
!      r_xy = sqrt((x(i) - x0) * (x(i) - x0) + (y(j) - y0) * (y(j) - y0))

!      if (r .eq. 0.0) r = eps
!      if (r_xy .eq. 0.0) r_xy = eps
!      sin_theta = r_xy / r
!      cos_theta = (z(k) - z0) / r
!      sin_phi = (y(j) - y0) / r_xy
!      cos_phi = (x(i) - x0) / r_xy

!      w = sqrt(r ** 2 + r0 ** 2 - 2.0 * r * r0 * sin_theta)
!      a_phi = 0.5 * q * a ** 2 * bt * exp(- w ** 2 / a ** 2)
!      br = a_phi * 4.0 * r0 * cos_theta / a ** 2 / r / sin_theta
!      btheta = a_phi * 2.0 * (r - r0 * sin_theta) / a ** 2 / sin_theta / r
!      bphi = a * bt / r / sin_theta * exp(- w ** 2 / a ** 2)

!!     rotated coordination y <-> z
!      bx(i, j, k) = bx(i, j, k) + sin_theta * cos_phi * br + cos_theta * cos_phi * btheta - sin_phi * bphi
!      by(i, j, k) = by(i, j, k) + sin_theta * sin_phi * br + cos_theta * sin_phi * btheta + cos_phi * bphi
!      bz(i, j, k) = bz(i, j, k) + cos_theta * br - sin_theta * btheta

      if (r .le. r0) then
         bx(i, j, k) = bx(i, j, k) - bt * e_dir * (y(j) - y0) / (r ** 2 + eps) * (r0 - r)
         by(i, j, k) = by(i, j, k) + bt * e_dir * (x(i) - x0) / (r ** 2 + eps) * (r0 - r)
         bz(i, j, k) = 0.0d0
      endif

      en(i, j, k) = pr / gmm1 + (bx(i, j, k) * bx(i, j, k) + by(i, j, k) * by(i, j, k) + bz(i, j, k) * bz(i, j, k)) / 2.0d0 +      &
         (mx(i, j, k) * mx(i, j, k) + my(i, j, k) * my(i, j, k) + mz(i, j, k) * mz(i, j, k)) / ro(i, j, k) / 2.0d0

   enddo
   enddo
   enddo

!----------------------------------------------------------------------------------------------------------------------------------|
   return
end subroutine bnd_user
